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We test the accuracy of our recently proposed empirical formula to model the recoil velocity im- 
parted to the merger remnant of spinning, unequal-mass black-hole binaries. We study three families 
of black-hole binary configurations, all with mass ratio q=3/8 (to maximize the unequal-mass con- 
tribution to the kick) and spins aligned (or counter aligned) with the orbital angular momentum, 
two with spin configurations chosen to minimize the spin-induced tangential and radial accelerations 
of the trajectories respectively, and a third family where the trajectories are significantly altered by 
spin-orbit coupling. We find good agreement between the measured and predicted recoil velocities 
for the first two families, and reasonable agreement for the third. We also re-examine our original 
generic binary configuration that led to the discovery of extremely large spin-driven recoil veloci- 
ties and inspired our empirical formula, and find reasonable agreement between the predicted and 
measured recoil speeds. 
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I. INTRODUCTION 

Thanks to recent breakthroughs in the full non-linear 
numerical evolution of black-hole-binary spacetimes [11, 
[3, Q , it is now possible to accurately simulate the merger 
process and examine its effects in this highly non-linear 
regime H S i, 0, B i E, [H E El M M, ll6|, [iTj, ll8| . 
Black-hole binaries will radiate between 2% and 8% of 
their total mass and up to 40% of their angular momenta, 
depending on the magnitude and direction of the spin 
components, during the merger [1, In addition, the 
radiation of net linear momentum by a black-hole binary 
leads to the recoil of the final remnant hole ^T^, '20', '21], 
M pi S Hi m, m, H i3 1 , 32, 33, M, 
40, l4li . which can have astrophysically important 

effects m Ha m, n, El, SI] . 

Merging black-hole binaries will radiate net linear mo- 
mentum if the two black holes are not symmetric. This 
asymmetry can be due to unequal masses, unequal spins, 
or a combination of the two. A non-spinning black-hole 
binary will thus only radiate net linear momentum if the 
component masses are not equal. However, the maximum 
recoil in this case (which occurs when the mass ratio is 
q « 0.36) is a relatively small - 175 kms'^ [H]. The 
complementary case, where the black holes have equal 
masses but unequal spins was first reported in [27[ and 
[2^. In the former case the authors calculated the re- 
coil velocity for equal-mass, quasi-circular binaries with 
equal-amplitude, anti-parallel spins aligned with the or- 
bital angular momentum direction, while in the latter 
case the authors used the same general configuration but 
varied the amplitude of one of the spins. In both the 
above cases the authors extrapolated a maximum pos- 
sible recoil (which is tangent to the orbital plane) of 
^ 460 kms~^ when the two holes have maximal spin. 
At the same time, our group released a paper on the 



first simulation of a generic black-hole binaries with un- 
equal masses and spins, where the spins were not aligned 
with the orbital angular momentum [2^. That configu- 
ration had a mass ratio of 1:2, with the larger black hole 
having spin a/m — 0.885 pointing 45° below the orbital 
plane and the smaller hole having negligible spin. The 
black holes displayed spin precession and spin flips and 
a measured recoil velocity of 475 kms~^, mostly along 
the orbital angular momentum direction. It was thus 
found that the recoil normal to the orbital plane (due 
to spin components lying in the orbital plane) can be 
larger than the in-plane recoil originating from either the 
unequal-masses or the spin components normal to the 
orbital plane. The maximum possible recoil arises from 
equal-mass, maximally spinning holes with spins in the 
orbital plane and counter-aligned. This maximum re- 
coil, which will be normal to the orbital plane, is nearly 
4000 kms"^ 

In (2^ we introduced the following heuristic model for 
the gravitational recoil of a merging binary. 

V^rccoii(g, di) = v,n ci + i;_l(cos(0 61 -f sin(^) 62) + v\\ e^, 

(1) 

where 



Vrr,. = A— \l+ B 



H- 



(i + g)V' 



(1 + 9)' 
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(2b) 



q 

?;|| = /^cos(e - 60)-— — ~ qd^l , (2c) 

1.2 X W kms-i [11], B = -0.93 here we 

flnd H = (6.9 ± 0.5) x 10^ kms-^, a, = SJm^, Si 
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and rrii are the spin and mass of hole i, q = mi/ 1712 is 
the mass ratio of the smaller to larger mass hole, the 
index _L and || refer to perpendicular and parallel to 
the orbital angular momentum respectively at the effec- 
tive moment of the maximum generation of the recoil 
(around merger time), 61,62 are orthogonal unit vectors 
in the orbital plane, and ^ measures the angle between 
the "unequal mass" and "spin" contributions to the re- 
coil velocity in the orbital plane. The angle Q was de- 
fined as the angle between the in-plane component of 
A = (mi +m2){S2/m2 — Si/mi) and the infall directjon 
at merger. The form of Eq. ([2a|) was proposed in j23i . [46ll . 
while the form of Eqs. (|2bp and ([2c|) was proposed in [23| 
based on the post-Newtonian expressions in In 
Rcf [48] we determined that K = (6.0±0.1) x 10'* kms'^ 
Although ^ may in general depend strongly on the con- 
figuration, the results of [13] and post-Newtonian calcu- 
lations show that ^ is 90° for headon collisions, and the 
results presented here indicate that ^ ~ 145° for a wide 
range of quasi-circular configurations. A simplified ver- 
sion of Eq. ^ that models the magnitude of Kccoii was 
independently proposed in 32], and a simplified form of 
Eq. m) for the equal-mass aligned spin case was proposed 
in ^9]. 

Our heuristic formula ([T]) describing the recoil veloc- 
ity of a black-hole binary remnant as a function of the 
parameters of the individual holes has been theoretically 
verified in several ways. In [4^ the cos0 dependence 
was established and was confirmed in [371] for binaries 
with larger initial separations. In Ref. [36| the decompo- 
sition into spin components perpendicular and parallel 
to the orbital plane was verified, and in [41'] it was found 
that the quadratic-in-spin corrections to the in-plane re- 
coil velocity are less than 20 kms^^. 

Consistent and independent recoil velocity calculations 
have also been obtained for equal-mass binaries with 
spinning black holes that have spins aligned/ counter- 
aligned with the orbital angular momentum [13, [1^ . Re- 
coils from the merger of non-precessing unequal mass 
black-hole binaries have been modeled in [3^. 

The net in-plane remnant recoil velocity arises both 
from the asymmetry due to unequal masses, which given 
its z — > — z symmetric behavior, only contributes to recoil 
along the orbital plane, and the asymmetry produced by 
the black-hole spin component perpendicular to the or- 
bital plane. Even if we can parametrize the contribution 
of each of these two components of the recoil in terms of 
only one angle, ^, the modeling of it appears in princi- 
ple very complicated. ^ may depend on the mass ratio 
(q) of the holes, as well as their individual spins and 
S2, but also on their orbital parameters such as initial 
coordinates and momenta, or initial separation and ec- 
centricity. We clearly have to reduce the dimensionality 
of this parameter space as part of the modeling process. 
In order to do so, we shall choose a model for ^ that 
only depends on q and for quasi-circular orbits. We 
then perform simulations to determine how accurately 
this reduced-parameter-space model for ^ reproduces the 



observed recoil velocities and find that ^ « 145°, inde- 
pendent of either g or A^ . 

The paper is organized as follows, in Sec. [IT] we re- 
view the numerical techniques used for the evolution of 
the black-hole binaries and the analysis of the physical 
quantities extracted at their horizons. In Sec. IIIII we re- 
view the post-Newtonian dynamics of binary systems in 
order to motivate our study of equivalent trajectories for 
unequal mass, nonspinning and spinning holes. We focus 
on four families of such configurations. In Sec. IIVI we 
give the initial data parameters for these families. The 
results of the evolution of those configurations are given 
in Sec. |Vl where we also introduce a novel analysis of the 
trajectories of the punctures and of the waveform phase 
to model the angle ^ in our heuristic formula Eq. ([1]). 
In Sec. IVII we analyze the generic configuration that led 
us to discover the large recoil velocities produced by the 
spin projection on the orbital plane of the binary. Here 
we use more refined tools to analyze the individual hole 
spins and momenta near merger time, when most of the 
recoil is generated. We end the paper with a Discussion 
section pointing out the need for further runs with higher 
accuracy to improve our first results, and an Appendix 
including the post-Newtonian analysis of the maximum 
recoil configuration. 



II. TECHNIQUES 

We use the puncture approach [49l] along with the 
TwoPuNCTURES [1^ thorn to compute initial data. In 
this approach the 3-metric on the initial slice has the form 
7ab = {fpBL + u)'^6ab, whcrc ipBL is the Brill-Lindquist 
conformal factor. Sat is the Euclidean metric, and u is 
(at least) on the punctures. The Brill-Lindquist con- 
formal factor is given by i/'sl = 1 + X]"=i '^i'/(2l^~ 
where n is the total number of 'punctures', is the 
mass parameter of puncture i (m^ is not the horizon 
mass associated with puncture i), and fi is the coor- 
dinate location of puncture i. We evolve these black- 
hole-binary data-sets using the LazEv_^1| implementa- 
tion of the moving puncture approach |3] . In our version 
of the mo ving puncture approach [2, |3] we replace the 
BSSN [5^. l53l. Im] conformal exponent (j), which has log- 
arithmic singularities at the punctures, with the initially 
field X = exp(— 40). This new variable, along with 
the other BSSN variables, will remain finite provided that 
one uses a suitable choice for the gauge. An alternative 
approach uses standard finite differencing of (j) [Sj] . 

We use the Carpet [55 , ,56] mesh refinement driver to 
provide a 'moving boxes' style mesh refinement. In this 
approach refined grids of fixed size are arranged about the 
coordinate centers of both holes. The Carpet code then 
moves these fine grids about the computational domain 
by following the trajectories of the two black holes. 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified 1-1- log lapse and a modified Gamma-driver 
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shift condition 0, [131, and an initial lapse a{t = 0) = 
2/(1 + The lapse and shift are evolved with 

dt-P'd^)a = -2aK (3a) 
dtP" = (3b) 
dtB" = 3/45tf''-r/B". (3c) 

These gauge conditions require careful treatment of x, 
the inverse of the three-metric conformal factor, near the 
puncture in order for the system to remain stable [3, 0, 
[13. In Ref. m it was shown that this choice of gauge 
leads to a strongly hyperbolic evolution system provided 
that the shift does not become too large. 

We use AHFiNDERDiRECT [sll to locate apparent hori- 
zons. We measure the magnitude of the horizo n sp in us- 
ing the Isolated Horizon algorithm detailed in [60] ■ This 
algorithm is based on finding an approximate rotational 
Killing vector (i.e. an approximate rotational symmetry) 
on the horizon, and given this approximate Killing vector 
the spin magnitude is 

{v^R'Kab)d^V (4) 

°^ J AH 

where Kab is the extrinsic curvature of the 3D-slice, (PV 
is the natural volume element intrinsic to the horizon, 
and is the outward pointing unit vector normal to 
the horizon on the 3D-slice. We measure the direction of 
the spin by finding the coordinate line joining the poles 
of this Killing vector field using the technique introduced 
in 3. Our algorithm for finding the poles of the Killing 
vector field has an accuracy of ~ 2° (see for details). 

We also use an alternative quasi-local measurement of 
the spin and linear momentum of the individual black 
holes in the binary that is based on the coordinate ro- 
tation and translation vectors [s^ . In this approach the 
spin components of the horizon are given by 

= ^U^'Kabd'V, (5) 

O"" J AH 

where — Sijdmkr™^^^'' , and r™ = a;™ — is the 
coordinate displacement from the centroid of the hole, 
while the linear momentum is given by 

= ^ / - K^ab)d^V, (6) 

OTT J AH 

where = ^e- 

We measure radiated energy, linear momentum, and 
angular momentum, in terms of V'4, using the formulae 
provided in Refs. [61I, [62|. However, rather than using 
the full Tpi we decompose it into i and m modes and 
solve for the radiated linear momenturn, d ropp ing terms 
with i > 5. The formulae in Refs. [6l|, [GJ] are valid 
at r = 00. We obtain highly accurate values for these 
quantities by solving for them on spheres of finite radius 
(typically r/M = 25,30,35,40), fitting the results to a 
polynomial dependence in I — l/r, and extrapolating to 



I = 0. We perform fits based on a linear and quadratic 
dependence on I, and take the final values to be the aver- 
age of these two extrapolations with the differences being 
the extrapolation error. 

We obtain a new determination of H in Eq. (j2bp 
using results from simulations performed by the 
NASA/GSFC m, PSU [H, and AEI/LSU [il^ groups. 
The simulations performed by these groups include runs 
with q ~ 1, and thus provide an accurate measure- 
ment of v±. We calculate H for each simulation (via 

H = vj_{al — af )(1 -I- q)^ /q^) and take the weighted av- 
erage {H) ± 5{H), where 

i 

6{X) = ^{X^) - (7) 

X is the quantity to be averaged, n is some specified 
power, and dXi is the uncertainty in a particular mea- 
surement of X. Note that we weight H and by the 
same w,. We find (H) = (6895 ± 513) kms^^. Fig- 
ure [1] shows the values of H obtained from each simula- 
tion as well as the average value of H. We can see that 
based on the AEI/LSU data, which take into account 
the initial recoil at the beginning of the full numerical 
simulations, one could fit linear corrections to H . How- 
ever, the deviations from H = const are only significant 

near D — q^ /{\ + q)^{a\ — qa\) = 0, when the spin- 
induced recoil is small (and hence the relative error in 
the spin- induced recoil is large). The absolute differences 
between the predicted and measured recoil velocities for 
the AEI/LSU results are within 20 kms"^ when we take 
H = 6895 kms-i. 



III. POST-NEWTONIAN ANALYSIS 

In order to compare results from the recoil due to un- 
equal masses and those due to spin effects as well, we will 
study systems with similar orbital trajectories. Since the 
radiated momentum due to unequal masses is a function 
of the orbital acceleration, these systems will all exhibit 
very similar unequal-mass contributions to the net recoil, 
which allows us to isolate the spin-induced contributions 
to the recoil. To generate families of binaries with simi- 
lar trajectories we use the formulae for the leading order 
post-Newtonian accelerations and choose configurations 
that minimize the effects of the spins on the trajectories, 
but have non-negligible spin contributions to the net re- 
coil. 

The relative one-body accelerations can be written 
as W\ 

d^dN + apN + d2PN + d,RR + dso + d,ss, (8) 

with 
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FIG. 1: The value of H calculated by inverting Eq. \2h\ 
as determined from simulations by the AEI, PSU, and 
NASA/GSFC groups. The thick line is the weighted aver- 
age and the thin lines are the expected uncertainty in the 
average. 
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(9c) 



aRR = < T"^ 



T 2 TO o 

18^2 + 25f2 

3 r 



6v^ - 2 15r' 
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(9d) 



aso = — J 6n[(n x v)i2S + —A)] - p x (75 + 3— A)] + 3f[n x (35 + —A)] 

7"' TO TO TO 



(9e) 



ass = "('^1 ' '^2) + 5i(n • ^2) + S2{n ■ Si) - 5h{n ■ Si){n ■82)}, 



(9f) 



where x = xi ~ X2, v — dx/dt, n = x/r, m — toi + TO2, /i — toiTO2/to, j/ = /i/to, (5to = toi — TO2, 5 = 5i + ^2, 
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and A = m(Sil'm2 — Si /mi), and an overdot denotes 
d/dt. 

The first four terms in Eq. ([8]) correspond to the 
Newtonian, first-post-Newtonian (IPN), second-post- 
Newtonian, and radiation reaction contributions to the 
equations of motion. The last two terms in Eq. ([H]) are 
the spin-orbit {SO) and spin-spin [SS) contributions to 
the acceleration. 

The radiated linear momentum due to the motion of 
the two holes has the form [471 



P 



N 



8 5m 
105 TO 



55^2 - 45r2 -I- 12- 



■ V 



TO 
I — 

r 



(10) 



plus higher post-Newtonian terms [63| , while the radiated 
linear momentum due to spin-orbit effects has the form 



P 



so 



8 ^^m 
15 



|4r(i/ X A) - 2v^{fL X A) 



{nxv) 3f{n-A) + 2{v-A) |. (11) 



Note also that the spin-spin coupling does not contribute 
to the radiated linear momentum to this order. 

In order to best study and model how the final remnant 
recoil velocity depends on the mass ratio and spins, we 
will chose configurations with black holes spinning along 
the orbital angular momentum. In this way the orbital 
plane will not precess and we can write [47| 



and 



S — Si -\- S2 — S^ 



rfi -|- ruX, 



(12) 



(13) 



where L^v = /^(a' x v) is the Newtonian orbital angular 
momentum, A = Lff x fi with Ljv = Ljv/|iAr|, and u = 
d(j)/dt is defined as the orbital angular velocity. 

Taking into account that the velocity remains in the 
orbital plane, i.e. Eq. psp . we find that the spin-orbit 
acceleration (Pe|) is given by 



1 



"SO - ^ 



55" 



3 — A-' 

TO 



2rS'X 



(14) 



and the radiated linear momentum is given by 
16 ^^toA" 



P 



so 



15 



-|(f2 „^2^2^^_ 2rrwn|, (15) 



and 
Pn 



8 Sm 2 (^\'^ 
105 TO \ r / 



hr^u^ - 2r^ +4— 
r . 



— ruj\ 



50r^w^ + 12r-' + 8— 
r . 



(16) 



Note that if we take the scalar product of these two in- 
stantaneous radiated momenta we obtain 

Pn ■ P'soKPnP^o) = ^os{(pn) 

= -^f{r, r, uj)/yj\g{r, f, uj){ll) 

where / = 4rf^ + (8to + 2Ar'^uj'^)P - r^uj^{4m + 25r^LJ^) 
and g{r,r,uj) — Ar^ — (16to — 12Ar^uj'^ — r)r*/r + 
{16m'^ + 2S2r^uj^m+1225r^uj'^ + 2r'^uj^)P/r^+uj^{6Am'^ + 
SOOr^uj^m + 2500r^uj'^ + r^uj"^). The fact that the factor 
of A" drops out of Eq. ((T7)) suggests that ^ (which is the 
angle between the cumulative radiated linear momenta) 
will depend only weakly on the spins through their affects 
on the orbital motion. Binaries with similar orbital tra- 
jectories should therefore have similar values for ^. Note 
that ^ may still be a strong function of trajectory and q. 

We will now turn to the question of identifying a subset 
of physical parameters of the binary that produce similar 
trajectories for unequal-mass, non-spinning and unequal- 
mass, spinning binaries in order to compare their recoil 
velocities and extract the spin contribution to the total 
recoil. 



A. similar radial trajectories 

An analysis of how ^ depends on configuration is 
greatly simplified if the trajectories of the spinning bina- 
ries are similar to the trajectory for a non-spinning binary 
with the same mass ratio. In order to accomplish this, 
we use the post-Newtonian expression for the spin-orbit 
induced acceleration Eq. p4)) . and choose configuration 
that minimize its effect. 

The radial component of the spin-orbit induced accel- 
eration wiU vanish ii 55" + 3^A" = 0. This leads to 

m 

the condition 

F= (3g + 2) + (3 + 2g)5 = 0, (18) 

where a = qai/a2 can take any positive or negative 
value. However, if we consider the algebraic average over 
the range < g < 1 at fixed F we find 



1 



{a) = -[a(g = 0)+a(g= 1)] = 



(19) 



5 

15' 6' 

and that a — (a) when q — S/8 (independent of F). 

We will thus study configurations with this mass ratio 
(which also produces a nearly maximum recoil velocity of 
« 175 kms""'^ for non-spinning unequal mass black hole 
binaries j23j). 

Hence the first family of black-hole-binary configura- 
tions that we will study is given by the choice 



F = 0, g = 3/8. 



(20) 



thus 



a2/ai - -g(3 + 2q)/{2 + 3q) = -9/20. (21) 

The total spin of the binary will in general be non- 
vanishing with 



57to2 = (a2 + q^ai)/{l + q)^ = 4a2/ll. 



(22) 
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TABLE I: Initial data parameters for quasi-circular orbits with orbital frequency u/M — 0.05. All sets have mass ratio 
q = m" /m" = 3/8. The 'F' series has a = 02/01 = -9/20 (hence F = qaxlai{2q + 3) + 3g + 2 = 0), and the 'S' series has 
S = 5*1 + S'2 = 0. The punctures are located along the x-axis with momenta Pi = (0, P, 0) and P2 = (0, — P, 0), and spins Si 
along the z-axis. m*' are the puncture masses, ra^ are the horizon masses. 



Config 


Q38 


F+0.2 


F_0.2 


F+0.4 


F-0.4 


S+0.64 


S-0.64 


XX IM 


-4.7455652 


-4.6889329 


-4.8008847 


-4.6310312 


-4.8548401 


-4.5310235 


-4.9561392 


X2/M 
Sl/M' 

siiW 


1.7604572 


1.8161037 


1.7042740 


1.8711650 


1.6475993 


1.897592 


1.6168224 


0.0000000 


0.015219622 


-0.015222140 


0.030437161 


-0.030447242 


0.048726127 


-0.048689700 


0.0000000 


-0.048702791 


0.048710847 


-0.097398914 


0.097431175 


-0.048726127 


0.048689700 


P/M 

J/M'^ 




n 1 n7n7Q9Q 


U.IUDOD / ^ ( 


n 1 07*^/19/1/1 


n 1 C\f\^A 7Q9 


u.iuo ( Do^y 


n 1 nf^Q9Qpisi 
u. luoyzyoo 


0.69498063 


0.6965546816 


0.6932382744 


0.6979616178 


0.6913258658 


0.6863415524 


0.7028440196 


0.69498063 


0.6630715129 


0.7267269819 


0.6309998643 


0.7583097987 


0.6863415524 


0.7028440196 


m\IM 


0.257487827988 


0.25319314 


0.253279647 


0.239665153 


0.239816706 


0.205915971 


0.206131153 


ml/M 


0.718534207968 


0.71621170 


0.716211394 


0.709030409 


0.70903488 


0.715832591 


0.715746409 


mx /M 


0.27582974 


0.27577886 


0.275791869 


0.275757065 


0.27577578 


0.2756959 


0.27558121 


/M 


0.73541100 


0.73541402 


0.735444371 


0.735334919 


0.735402505 


0.7351861 


0.734888095 


al 


0.000 


0.20012582 


-0.20013825 


0.4002982516 


-0.400383662 


0.6411766 


-0.64119084 


a\ 


0.000 


-0.090053523 


0.0900619119 


-0.180130779 


0.18015874 


-0.090153 


0.09015883 




1.00001 


1.00001 


0.999997 


1.00001 


0.999997 


1.00001 


0.999991 



B. similar tangential trajectories 

We can also choose a configuration where the tangen- 
tial component of the acceleration due to the spin-orbit 
coupHng vanishes, i.e. 

= SI + S'^=Q. (23) 

This translates into the condition 

a2/ai = -q^ = -9/64 (24) 

when q = 3/8. Note that now, the radial acceleration, as 
parametrized by F, is non vanishing 

F = (3q + 2) + (3 + 2q)qai/a2 - -55/8. (25) 

Thus, for (7 7^ 1, we cannot make both the radial and tan- 
gential components of the spin-orbit acceleration vanish 
at the same time by a simple choice of physical parame- 
ters of the binary. 



eters that we will denote by Q, F, S, and A. The Q-series 
has initially non-spinning holes, the F-series has F ^ Q 
(See Eq. (|18p ): hence zero PN spin-orbit-induced radial 
acceleration, the S-series has total spin 5 = 0; hence 
zero PN spin-orbit-induced tangential acceleration, and 
the A-series has neither F = Q nor 5* = 0; hence both 
spin-obit-induced accelerations are non-vanishing. The 
puncture masses were fixed by requiring that the total 
ADM mass of the system be 1 and that the mass ratio of 
the horizon masses be 3/8. The initial data parameters 
for these configurations are given in Tables U and [TTl We 
obtained initial data parameters by choosing spin and 
linear momenta consistent with 3PN quasi-circular or- 
bits for binaries with mass ratio q — i/S and then solve 
for the Bowen-York ansatz for the initial 3-metric and 
extrinsic curvature. This method was pioneered by the 
Lazarus project (63| (See Fig. 35 there), and then used 
in the rest of the breakthrough papers [S, 7, 8, 39, 48, 6^ 
by the authors (in Ref. [11] we used the PN expressions 
for the radial component of the momentum as well). 



IV. INITIAL CONFIGURATIONS 



We choose quasi-circular initial configurations with 
mass ratio q — mi/m2 = 3/8 from four families of param- 



V. RESULTS 

We evolved all configurations given in Tables U and |TT] 
using 10 levels of refinement with a finest resolution of 
h — M/80 and outer boundaries at 320M except configu- 
ration A-)_o.g, where we used an additional coarse level to 



push the outer boundaries to 640M. In all cases, except 
where noted otherwise, we set the free Gamma-driver pa- 
rameter in Eq. ^ to ry = 6/M 

In a generic simulation both the unequal mass and spin 
components of the recoil are functions of the trajectory. 
To single out each individual effect we perform runs cho- 
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TABLE II: Initial data parameters for quasi-circular orbits 
with orbital frequency LojM = 0.05. All sets have mass ratio 
q = /m2 ~ 3/8. The punctures are located along the 
X-axis with momenta Pi = (0, P, 0) and P2 = (0, — P, 0), and 
spins Si along the z-axis. are the puncture masses, are 
the horizon masses. In this series neither F nor S vanishes. 



Config 


A-j-0.9 


A_o.9 


xi/M 


-4.5443438 


-4.8662563 


X2/M 
Sl/M^ 


1.573114 


1.9275192 


0.0000000 


0.0000000 


0.48873779 


-0.48581609 


P/M 

UjM^ 

J/M'^ 


0.10276465 


0.11089309 


0.6286584770 


0.7533827395 


1.117396265 


0.2675666537 


m\/M 


0.2545666 


0.2545806 


ml/M 


0.2822299 


0.284150275 


mf/M 


0.2733564 


0.2726292 


/M 


0.728824 


0.7270093 


a\ 


0.0000000 


0.00000 


a| 


0.920196524 


-0.9192121 


Madm/M 


1.000000 


0.999991 



sen to follow similar trajectories. In order to compare 
recoil velocity directions between these runs we need to 
rotate each system so that the final plunge (where most of 
the recoil is generated) occurs along the same direction. 
We do this in two ways. First, as demonstrated in Fig.[2l 
we plot the puncture trajectory difference r ~ — r2 
(where ^(t) is the coordinate location of puncture i at 
time for each case and rotate the trajectories by an 
angle <I>tiack so that they all line up with the Q38 tra- 
jectory during the late inspiral and merger phases. Note 
that by taking the differences between trajectories we re- 
move effects due to the wobble motion of the center of 
mass. Second, we measure the phase of the dominant 
(£ = 2, m = 2) mode of -04 at the point of peak ampli- 
tude and take half the phase difference between each case 
and Q38 (a rotation of about the z-axis will introduce 
a phase difference of —20 in the m = 2 components of 
1/^4) . We denote this latter rotation angle by $^4. We 
get reasonable agreement between these two measures of 
the rotation angle (See Tabic |TTT|. This type of rotation 
may also be needed when comparing results from differ- 
ent resolutions of the same configuration (i.e. when the 
phase error, but not the amplitude error, is large). In 
Table IIVI we give the components of the recoil velocity 
for a set of Qsg simulations with t] = 2/M. This value 
of 77 leads to a poorer effective resolution than for our 
standard choice of 77 = 6/M. Consequently there is a 
relatively large phase error in the low resolution results. 
After performing the rotation, the recoil velocities agree 
to within errors. 

Note that there is no rotation which will make the 
or A_o.9 trajectories line up with the Qss trajec- 
tory. In these cases the hangup-effect [1] due to spin-obit 
coupling significantly alters the trajectories (See Fig. [3]). 

Once we have found the correct rotation angle we ob- 



TABLE III: The rotation angle needed to align the trajecto- 
ries of each simulation with the Qag simulation as measured 
directly from the orbital trajectories ("I>track) and using the 
phase of the waveform at the point of maximum amplitude 
($^4). Note that provides the rotation angle modulo 
180°. 



Config 


'I'track 


*V.4 


^track - 


F+0.2 


25° 


34.5° 


9.5° 


F_0.2 


-28° 


-35.5° 


7.5" 


F+0.4 


56° 


63.1° 


7.1° 


F-0.4 


-44° 


-40.0° 


4.0° 


S+0.64 


5° 


9.7° 


4.7° 


S_0.64 


56° 


44.6° 


11.4° 


A+0.9 


* * * 


12.3° 


* * * 


A-0.9 


* * * 


-15.9° 


* * * 



^4 



TABLE IV: Results of the recoil velocity for the Qag con- 
figuration with 77 = 2/M at two different resolutions. After 
correcting for the phase error, equivalent to a rotation, the 
two recoils agree. Here 'Rtrack' denotes the velocity after ro- 
tating by the angle $track and 'P^^' denotes the velocity after 
rotating by the angle . 



h 


'I'track 




K 


Vy 


M/80 


34° 


36.5° 


-163 ± 12 


-46± 11 


M/80 Pv,4 


+ + + 


+ + + 


-103 ±12 


-134 ± 11 


M/80 Rtrack 


+ + + 


+ + + 


-109 ± 12 


-129 ± 11 


M/lOO 








-109 ± 14 


-133 ± 12 



tain ^ via 



i?[$]^rccoil, 

1 ~ Vq38, 
COS(0 = Kipin • Vq38, 



^pin — 



(26) 



where t^ecoii is the measured recoil velocity, rotates 
Kocoii by an angle $ in the xy plane, and Vgag is the recoil 

of the Q38 configuration. Note that when a\ — qaf < 
we need to replace ^ by tt — ^ in formula since the 
coefficient w_l in Eq. ^ is negative. We calculate two 
different values of ^, Attack and ^^4, based on the rotation 
angles $tiack and respectively. We obtain an addi- 
tional measurement of ^ by solving for cos ^ using Eq. ([T]) 
and the measured values of the recoil magnitude. We de- 
note this latter measurement of ^, which is unaffected by 
rotations, by ^Formula, where 



Formula 



cos 



i{qY - w_L(g, aill,a2ll) 



Ih2 



2wm('7) u_l((7, ai'l,a2ll) 



(27) 

Vyn is given by Eq. ([2a|) , v± is given by Eq. (j2bp , and v is 
the measured magnitude of the recoil velocity. 

We summarize the results of our simulations in Ta- 
bles |V] and IVII All configuration, with the exception 
of the 'A' series, have radiated energies in the range 
Erad/M = 0.021 ± 0.002 and radiated angular momenta 
in the range Ji-ad/M'^ = 0.15 ± 0.01, which is consistent 
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FIG. 2: The trajectory differences r = ri — r2 for the 'F' and 'S' series rotated so that the late-inspiral matches the Q38 
trajectory. The plots show the rotation angle $track 
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with these trajectories being essentially the same for all configurations (See Fig. [2). 
I 

I 



We obtain weighted averages for ^ for the 'F' and 'S' series of (^track) = (152 ± 9)°, (^^,) = (143 ± 14)°, and 
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TABLE V: The recoil velocities (prior to any rotation), radiated energy and angular momentum, and ^ for the 'Q' and 'F' 
series, ^track is calculated using <l?track and Eq. (|26p . is calculated using "1>^^ and Eq. (|26p . ^Formula is calculated from 
the given recoil magnitude using Eq. ((27}. It^trackl' I^vm"'*!' ^^'^ l^avg""*! B.Te the recoil velocities as predicted by Eq. ^ with 
C = Ctrack, C = Ci^4, and ^ = (^) respectively. 



Config 


Q38 


F+0.2 


F-0.2 


F+0.4 


F_0.4 


Jrad/M^ 


0.0210 ± 0.0003 


0.0202 ± 0.0003 


0.0219 ±0.0004 


0.0193 ±0.0002 


0.0228 ± 0.0004 


0.1503 ± 0.0030 


0.1471 ±0.0005 


0.1576 ±0.0015 


0.1399 ±0.0016 


0.1625 ±0.0010 


y"[kms-i] 


-94± 11 


-177 ± 10 


-15± 14 


-223 ± 12 


15± 14 


y^[kms-^] 


-141 ± 5 


-85± 12 


-155 ±5 


33± 18 


-127 ±4 


\V\[kms-^] 


169.5 ± 7.4 


196.4 ± 10.4 


155.7 ±7.2 


225.4 ± 12.2 


127.9 ±4.3 


Ctrack[deg] 

5*4 [deg] 


0° 
0° 


(143 ±31)° 
(154 ±43)° 


(178 ± 73)° 
(127 ±41)° 


(147 ± 20)° 
(173 ±25)° 


(169 ±21)° 
(179 ±21)° 


CFormula[deg] 


0° 


(127 ± 26)° 


(131 ± 15)° 


(134 ±20)° 


(144 ±6)° 


IK?rkl[kms-i] 


175 


202 ±9 


142 ±3 


232 ± 10 


112 ±9 


|F^7<^|[kms-il 
|FP'g°''|[kms-i] 


175 


205 ±9 


158 ±21 


240 ± 5 


110 ±5 


175 


203 ±3 


150 ±4 


231 ± 5 


127 ±8 



TABLE VI: The recoil velocities (prior to any rotation), radiated energy and angular momentum, and ^ for the 'S' and 'A' 
series. Note that although we report the calculated values for ^ based on for the 'A' series, here ^ is not well defined 
because the unequal mass component of the recoil is not given by the Qas recoil, ^track is calculated using $track and Eq. (|26|l , 
is calculated using and Eq. p6|l . ^Formula is calculated from the given recoil magnitude using Eq. ([27|l . I^^rlckl' I'vl^'^li 
and |yavg°'*| are the recoil velocities as predicted by Eq. ([1} with ^ — ^track, C = C*4j and ^ — (^) respectively. 



Config 


So+0.64 


S-0.64 


A+0.9 


A_o.9 


Jrad/M^ 


0.0209±0.0003 


0.0203 ± 0.0002 


0.050668 ± 0.000974 


0.01274 ± 0.00003 


0.152 ±0.0007 


0.146 ±0.001 


0.2999857 ± 0.00708 


0.092 ±0.001 


V"[kms-i] 


-122 ± 18 


-119±5 


13 ±30 


48 ±24 




-181 ± 15 


31 ±4 


-63 ±2 


-340 ± 8 


ll^Kkms"'] 


218.3 ± 16.0 


123.0 ±4.9 


64.1 ± 5.9 


343.4 ± 8.6 


5track[dcg] 

5*4 [deg] 

5Formula[deg] 


(160 ±31)° 
(142 ±28)° 
(124 ±22)° 


(148 ± 11)° 
(137 ±7)° 
(150 ±7)° 


* * * 

(158 ±7)° 
(159 ±2)° 


* * * 

(93 ± 7)° 
(149 ± 19)° 


IK?rckl[kms-i] 


237 ± 10 


125 ± 10 


* * * 


* * * 


IF^^'^likms-i] 
IKP-'^likms-i] 


230 ± 16 


135 ±7 


68 ±22 


259 ± 18 


231 ± 5 


127 ±8 


108 ± 28 


340 ±9 
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TABLE VII: The average value |(f)| of the trajectories for 
each configuration. The larger the value of |(f)| the slower 
the inspiral. 

Config \{r)\ Config Config \{r)\ 

038 0.85 8303 F+0.2 0.902234 F-q.2 0.7982157 

F+0.4 0.936172 F-oa 0.76745 5+0.64 0.845 197 

5-0.64 0.95333 A+o.g 0.365654 A-o.9 1.39869 



(^Formula) = (144 ± 7)°, where we use Eq. (O to obtain 
the weighted average and uncertainty. These weighted 
averages are consistent with the measured values of ^. 
The weighted average over all three measurements of ^ is 
{£_) — (145 ±10)°. Interestingly, (^) provides an accurate 
prediction for the recoil velocity of the A_q g configura- 
tion. This result is unexpected because the recoil due 
to unequal masses is a function of the mass ratio and 
the trajectories (i.e. the accelerations of the masses over 
time). For the 'F' and 'S' configuration the trajectories 
are very similar to Qss, and hence the unequal mass com- 
ponents of the recoil are expected to be very similar to 
Qss- However, the spin-orbit coupling induced hangup 
effect in both ^+0.9 and ^-0.9 greatly affects the trajec- 
tories (See Fig. (3]), as well as the radiated energy and 
angular momenta. If we consider the radiated linear mo- 
mentum averaged over an orbit, then we see that the 
slower the inspiral (i.e. the closer to a closed orbit), the 
smaller the average recoil. Hence we expect that ^+0.9 
will have a smaller unequal-mass recoil than Qssi while 
yl_o.9 will have a larger one. To quantify how much the 
orbits close we take the average of f = ri — r2 over the 
trajectory from the beginning of each simulation until 
|r| ~ 0.1. The resulting averages |(r)| for the 'Q', 'F', 
'S', and 'A' families are given in Table IVTIl The mean 
and standard deviation of \{r)\ for the 'Q', 'F', and 'S' 
configurations is \{r)\ = 0.865 ± 0.070. The A+0.9 and 
A-0.9 configuration lie 7.1cr and 7.6a below and above 
this mean respectively, while the results for the other 
configurations lie within lAa of the mean. 

As seen in Fig. U the angle ^ appears, at least quali- 
tatively, to be independent of A. This is in agreement 
with our post-Newtonian analysis in Eq. (flT)) . It is also 
consistent with our intuition that similar trajectories im- 
ply similar angles between the unequal-mass and spin 
contributions to the recoil, and it seems that the small 
differences in the trajectories produce some scatter on 
the values, but this is apparently mostly due to the nu- 
merical error generated during the simulations. It would 
be interesting to use this value of ^ to model the recoil 
velocity distribution in galaxies. 



VI. GENERIC EVOLUTION REANALYZED 

In light of our new understanding about the modeling 
of the recoil velocity, we re-examine our original generic 
binary configuration, which we denote by SPG. The SP6 
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FIG. 4: ^ versus A/m — S^jtn^ — Si/m\ as calculated in this 
work for a mass ratio 5 = 3/8 and from the data published 
by the NASA/GSFC group for a mass ratio q — 2/3 provided 
in Ref. fsl]. We plot ^tiack, C«/'4i and ^Formula for the 'F' and 
'S' configurations and ^Formula for 'A' configurations. The 
thick horizontal line and the two thin horizontal line show 
the average value (5) and its uncertainty (as calculated in 
this work from our simulations). The data displays significant 
scatter, but appears to be consistent with ^ = const. 



configuration has a mass ratio of g = 1/2 with the larger 
hole having specific-spin a/m — 0.885 and spin pointing 
45° below the orbital plane, and the smaller hole having 
negligible spin. We also evolved a similar configuration, 
which we will denote by SP6R, that is identical to SP6, 
but with the spin rotated by 90° about the z-axis. We 
evolved both configurations using the same grid structure 
as in the previous section, but used tj — 2/M rather than 
6/M. This choice of smaller rj has the effect of reduc- 
ing the effective resolution, but makes calculations of the 
quasi- local linear momentum and spin direction more ac- 
curate (See Ref. [1^) by reducing coordinate distortions. 
The initial data parameters for the two configurations 
are given in Table IVIIII The drop in effective resolution 
when reducing 77 from 6/M to 2/M is significant. In our 
simulations we found that a Q38, t] — 2/M run with cen- 
tral resolution of M/lOO had a slightly larger waveform 
phase error than an equivalent A//80 resolution run with 
7] = 6/M, while an M/80 run with ?/ = 2/M displayed 
a significant phase error. We have found in general that, 
with our choice of gauge, the coordinate dependent mea- 
surements, such as spin and linear momentum direction, 
become more accurate as 77 is reduced (and /i — > 0). How- 
ever, if 77 is too small {rj ^ 1/M), the runs may become 
unstable. Similarly, if 77 is too large {rj > 10/M), then 
grid stretching effects can cause the remnant horizon to 
continuously grow, eventually leading to an unacceptable 
loss in accuracy at late-times. We have found that a value 
of 77 = 6/Af provides both very high accuracy in the com- 
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TABLE VIII: Initial data parameters for the SP6 and SP6R 
configurations, nip is the puncture mass parameter of the two 
holes. SP6 has spins Si = (0, S, -5") and 5^2 = (0,0,0), mo- 
menta P = ±(Pr,P±,0), puncture positions xi = {x+,d,d) 
and X2 ~ {x-,d,d), masses mi and m2, and Maom/M — 
1.00000 ± 0.00001. SP6R has the same parameters as SP6 
with the exception that Si = ( — 5,0, —S). 



rUp/M 0.3185 
X+/M 2.68773 
X-/M -5.20295 



d/M 0.0012817 mi/M 0.6680 
Pr/M -0.0013947 mz/M 0.3355 
P±/M 0.10695 S/M'^ 0.27941 



puted waveform at modest resolutions, while keeping the 
remnant horizon size nearly fixed at late-times. 

We measure a net recoil of Kecoii = 375±18 km s^^ and 
K-ccoii = 848 ± 20 kms-i for SP6 and SP6R respectively. 

The analysis of the recoil in SP6 and SP6R is compli- 
cated by the fact that the orbital plane precesses signifi- 
cantly during the merger. Thus, we cannot associate the 
xy components of the recoil with the in-plane recoil (as 
was done tentatively in Ref. [13 )• order to measure 
the precession of the orbital plane we need an accurate 
measurement of the orbital angular momentum. Here we 
use the approximate formula 



^orbit 



X P, . 



(28) 



where fi is the coordinate location of puncture i and Pi is 
the quasi- local momentum 39] , given by Eq. ^ , of black 
hole i. In Fig. [5] we show the orbital angular momentum 
of the SP6 configuration versus time. Note the rapid 
change in direction near merger (a common horizon was 
first detected a.t t = 207. 4M), and as seen in Fig. [51 most 
of the recoil is generated about 3M to 30M after merger 
(here we assume that waveform features seen at t — t 
for an observer at r = 40Af were generated by dynamics 
near the horizons a.t t t ~ 40M). This rapid change 
in direction has a strong effect on the computed recoil 
due to the cos Q and cos^ dependence of Wrccoii- That is, 
rapid physical changes in the orbital plane and spin direc- 
tion, lead to relatively large errors in the direction (but 
not magnitude) of both the spin and orbital angular mo- 
menta when the resolution is below some threshold. This 
in turn, leads to relatively large errors in the measured 
recoil. Thus it is not surprising that this new calculation 
of the recoil velocity for SP6 is 100 kms~^ smaller than 
the value we reported i n I28l (note that we used a higher 
effective resolution in '28'|, thus we expect those values 
to be more accurate). These large errors will not be ob- 
served in more symmetric binaries where cither the spin 
or angular momentum axes are fixed. 

We can obtain an approximate measurement of a|| and 
a± using Eq. (^5)) and the measured direction of the 
spin. This estimation is only approximate due to the 
coordinate dependent nature of both calculations. We 
find that for SP6, a|| and a± vary little over the course 
of the run with values at merger of a» = —0.62 ± 0.03 
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FIG. 5: The normalized orbital angular momentum vector 
£ — L/\L\ versus time for the SP6 configuration up to merger. 
Note the rapid change in the direction at late times. 



and a± — 0.62 ± 0.03 (which are within errors of the 
initial values). However, the SP6R configuration does 
show a definite change in a over time, with merger val- 
ues of a|| = -0.69 ± 0.03 and aj. = 0.54 ± 0.03. We 
can use Eq. ^ to give estimates for the predicted re- 
coil velocity if we make the following assumptions: (1) 
^ = (^), (2) e for SP6R is rotated by tt/2 radians with 
respect to SP6, and (3) Qq is the same for SP6 and 
SP6R. Given these assumptions and the above range 
of the values for a|| and a±, we can perform a non- 
linear least-squares fit of the recoil magnitude for SP6 
and SP6R to obtain Oo. The resulting predictions for 
the recoil magnitude are Vspe — (500 ± 60) kms~^ and 
VsP6R — (1120 ± 130) kms~^. Both predictions are 
within 2a of the actual measured values and have an ab- 
solute error of 32%. If we fix ay and a± to their average 
values and vary our guess for ^ over the range (0, 360°), 
we find that the predicted values for Vsp6 a-nd Vsp6R lie 
in the ranges (462, 495)kms"i and (1048, 1120)kms-i re- 
spectively. 

The SP6 configuration demonstrated that the in-plane 
component of the spin can be the dominant contribution 
to the recoil. Given this observation, it becomes very 
important to accurately model this recoil. In Appendix El 
we derive a post-Newtonian model for the recoil produced 
by this in-plane component and show that it predicts the 
cos Q dependence in our empirical formula. 



VII. DISCUSSION 

Interestingly, most of the recoil velocity imparted to 
the remnant is generated at around merger time (more 
precisely, as seen in Fig. [SI within the first few tens of 
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enough accuracy for astrophysical applications. In par- 
ticular we have seen that the determination of an average 
value for the angle ^ of 145° seems to work not only for 
the F and S sequences, but also when we move off of 
these sequences towards more generic binaries. However, 
the formula should definitely be used with caution in an 
untested regime, especially when the trajectories are sig- 
nificantly altered by spin-orbit effects. 
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FIG. 6: The recoil speed (V = |l^|) for the SP6 configuration 
as measured from tpi a.t r = 40M as a function of time, as 
well as the time derivative of the recoil speed {dV/dt = V ■ 

V), and the magnitude of ip4,. Here the initial data burst is 
excluded from the calculation. Note that peak in dV/dt is 
located between t — 250M and t = 270M and occurs about 
2M latter than the peak in \^4\- A common horizon was 
first detected at t = 207. 4M, strongly suggesting that most 
of the recoil velocity is built up around merger time (since 
the observer is at r = 40M, features in the waveform at time 
t = T originated near the horizon(s) at time t ^ t — 40Af). 

M after merger. See also Refs. [H [H H^-), a non- 
linear regime where post-Newtonian approximations are 
not expected to work, but where the 'Lazarus' approach 
[H, ia H, 113, [ei. Hi] can be successfully applied 

Although an accurate modeling of ^ is challenging, 
starting from an ansatz that ^ = ^(g. A), we have found 
that, for quasi-circular orbits, ^ is qualitatively indepen- 
dent of either A or g for g = 3/8, g = 2/3 (based on the 
results of Ref. [Ill), and g = 1/2 (based on SP6). Note 
that the ^ that we measure is consistent with a simi- 
lar parameter introduced in Ref. [s^ . where they found 
^ = 147° (in our notation), based on a least-squares fit 
of the magnitude of the recoil versus a simplified version 
of Eq. ([T]). We know from the results for headon colli- 
sion (where ^ = 7r/2), that ^ is a function of eccentricity. 
However, for quasi-circular orbits, it appears to vary only 
marginally with either q or A. Further long-term simu- 
lations with high-accuracy (including extrapolations to 
/i — > and 77 — > 0) and further separated binaries will be 
needed in order to obtain a highly accurate model for ^. 
In particular, the — > limit will be important because 
the recoil depends sensitively on the linear momenta and 
spin directions of the individual black holes near merger 
(where gauge effects are most severe), and hence we need 
to take the 77 — > limit in order to accurately measure d?, 
L, and O. Nevertheless, our simple formula holds with 



APPENDIX A: POST-NEWTONIAN MODELING 

Here we provide a brief post-Newtonian analysis of 
the configurations that maximize the recoil velocity for 
spinning black holes. The spin-orbit-coupling (SO) con- 
tribution to the radiated linear momentum is given by 
Eq. HI]). 

We will restrict our analysis to planar orbits. Hence 
we have 

V = rii ~\- rwA, (-A-1) 

where A — Ln x h, Ljy = Ln/\Ln\, lu is the orbital an- 
gular velocity, and Lj^ = fi(x x v) is the Newtonian or- 
bital angular momentum. We shall take L^r = z. Hence 

\ ~ z X n and nx\ = nx{zxh)— z (A2) 

We observe that the third and fourth terms in Eq. (fTTj) 
only contribute to the recoil along the z-axis since 

n X V ~ rujz (-^-3) 

This contribution to the recoil velocity might well be the 
leading one, hence, in order to maximize the total recoil 
we seek to align, as much as possible, the first two terms 
in Eq. (jlip with the z-axis. This is achieved by having 
the spin of the black holes lie in the orbital plane, i.e. 

A = A„n + AaA. (A4) 

We then explicitly obtain the following products 

vxA = (t-Aa - rcjA„)z, (A5) 

n X A Axz, (A6) 

v-A = rAn+rLjAx. (A7) 
Plugging this into Eq. (fTTI) we find 

P^^ = -A^{(2r2 -4r2c.2)AA - 9rra;A„}z. (AS) 
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This clearly displays the fact that the recoil will 
be maximized when A takes the maximum magnitude 
(equal mass and opposite maximally rotating black holes) 
and varies sinusoidally with its projection along the line 
joining the holes. Note that if we define the angle be- 
tween h and A as we can write the above equation 
as 

i II 

= A(r) I A| cos 6* + B(r) I A| sin 61 

= C(r)|A|cos(0-0o(r)). (A9) 

This cos^ dependence in the recoil was the motivation 
for proposing the now- verified cos Q dependence in our 
empirical formula Eq. iPc]) for the recoil. 

Note that this analysis applies to the radiated linear 
momentum flux. Hence we have assumed that the larger 
the radiated linear momentum flux, the larger the total 
radiated linear momentum. 

It is also interesting to see if the unexpectedly large 
magnitude of the maximum out-of-plane recoil, compared 
to the in-plane recoil, can be understood using the post- 
Newtonian expression for the radiated linear momentum. 



i.e. Eqs. ([T5]) and (See Ref. [Ill for a similar anal- 

ysis). To do this, we used the post-Newtonian formulae 
for the radiated linear momentum along with the nu- 
merical trajectories for runs with the spins in the plane 
and perpendicular to the plane. We found that the post- 
Newtonian formulae predicted that the maximum out- 
of plane recoil will be approximately twice (almost 9/4) 
as large, rather than (the observed) « 8 times as large, 
as the maximum in-plane recoil. Thus we see that the 
magnitude of the out-of plane recoil arises from nonlin- 
ear dynamics at merger not fully captured by the post- 
Newtonian formalism. One may then conclude that, 
while the post-Newtonian approximation gives the cor- 
rect dependence of the recoil on the physical parameters, 
such as the scaling of the recoil velocities with the com- 
ponents of the spins parallel and perpendicular to the 
angular momentum, it is much less accurate when de- 
scribing the amplitude of the recoils. Thus we find that 
post-Newtonian formalisms provides the correct form for 
our semi-empirical formula ([1]), but does not provide ac- 
curate measurements of the magnitudes of the constants 
in that formula. 
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